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Abstract 

Pattern formation on semiconductor surfaces induced by low energetic ion-beam erosion under 
normal and oblique incidence is theoretically investigated using a continuum model in form of 
a stochastic, nonlocal, anisotropic Kuramoto-Sivashinsky equation. Depending on the size of the 
parameters this model exhibits hexagonally ordered dot, ripple, less regular and even rather smooth 
patterns. We investigate the transitional behavior between such states and suggest how transitions 
can be experimentally detected. 

PACS numbers: 68.55.-a,79.20.-m, 02.60.Lj 
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Self-organized structure formation on the nanoscale induced by ion-beam erosion, i.e. 
the removal of target material by bombarding its surface with ionized particles, has recently 
turned into a highly active research area of surface science with deep connections to the 
modern theory of non-equilibrium systems. Since the work of Navez et al. [1] in 1962 it is 
known that, under oblique incidence of the ions, often washboard-like ripple patterns ori- 
ented normally to the beam direction can be observed on the surface (for recent reviews 
cf. 

S3). 

For normal incidence of the ions, however, just a sort of rough, unstructured 
surface evolution had been expected. So, it came as a surprise when Facsko et al. 3] 
observed quite regularly ordered, hexagonally arranged dot structures for semiconducting 
GaSb under normal incidence of low energetic Ar + -ions. Parallel to that, Frost et al. 5( 
found for rotated InP, InAs, InSb and GaSb targets under oblique ion incidence a variety of 
distinct patterns such as dot or square structures and even rather flat or smooth surfaces. 
Rather recently, experiments under oblique incidence without target rotation have revealed 
an even more puzzling picture. Using Si targets and ion energies < 20006^, Ziberi et al. 
0] experimentally detected transitions from ripples oriented normally to the beam direction 
to rather smooth surfaces by increasing the incidence angle. Similar results have also been 
obtained by Ziberi et al. |7[ , for Ge and Si targets, where in the latter case almost perfectly 
straight ripples have been identified. Finally using GaSb targets and low ion fluxes, Allmers 
et al. ,8j identified a transition from hexagonal dot patterns to ripples oriented tangentially 
to the beam direction for small off-normal incidence. 

The theoretical understanding of the hexagonally ordered dot structures under normal inci- 
dence has posed a particular challenge because the standard continuum model for the evolu- 
tion of the morphology H (x, t) for ion-beam erosion, the anisotropic Kuramoto-Sivashinky 
equation (aKSE) 

ESQ does not seem to reproduce such patterns in its isotropic limit 
(iKSE), d t H = Fq + oxV 2 !! + a 2 V A H + a 3 (VH) 2 + 77. Therefore, two distinct generalizations 
of the iKSE have been put forward: Using different physical reasoning, Kim et al. 
and Castro et al. |3] have suggested the inclusion of a term proportional to V 2 (ViJ) 2 . Al- 
though this extended KSE (already studied in the context of amorphous surface growth Q|) 
might show patterns resembling hexagonally short-range ordered dots for some parameter 
values Q], it typ ically exhibits an irregular cellular pattern with remedies of a hexagonal 
arrangement [lj]. As an alternative to explain the hexagonally ordered dots, Facsko et al. 
jl^ suggested the inclusion of a damping term bH with b < in the iKSE that is known to 



2 



lead to hexagonal structures. This term, physically interpreted as the effect of redeposition 
of sputtered particles, fails to fulfill the fundamental symmetry of translation invariance 
in erosion direction H — > H + z with z =constant. To fix that, Facsko et al. suggested 
to replace bH by b(H — H) with H being the spatial average over some s amp le area. In- 
terpreting H as the spatial average over the whole sample, we showed in that such a 
non-local iKSE can be rigorously rewritten as a (local) damped iKSE (idKSE) by means of 
a temporally nonlocal transformation and, therefore, put the applicability of the idKSE on 
solid theoretical grounds. 

As a co_ of H Q, severa, im po rt a„ t Rations arise: (i) How does a m in imal 



model for ion-beam erosion under normal and oblique angle of incidence look like and what 
are its eminent emergent patterns? (ii) How does the transition from dots to ripples with 
varying angle of inclination happen? Is it a gradual change or a bifurcation? Does it occur at 
a zero or non-zero angle of incidence? (iii) How do experimentally observed smooth patterns 
with very low surface roughness fit in such a description? (iv) How do generic transitions 
between distinct patterns as function of the angle of incidence look like? In this letter, we 
theoretically investigate these questions on the basis of an anisotropically generalized version 
of a nonlocal KS equation and its damped counterpart. 

Model equation. - Generalizing the recipe in 16|, the balance equation for the evolution 
of the height function H(x,y,t), d t H = V ■ Jh + F + rj, is determined by (i) relaxational 
currents caused by surface diffusion and assumed to be dominantly isotropic modeled by 
V- Jh = a 2 V A H with a 2 < and (ii) and detachment contributions F = F[d x H, d y H, H — H] 
that depend on partial derivatives of H, H — H and all admissible combinations of them. 
Possible anisotropies stemming from ion-induced surface diffusion are at this stage neglected. 
Orienting the coordinate system such that the tangential component of the incoming ion 
flux is parallel to the x-axis, any term depending on odd order derivatives with respect to y 
is excluded due to the reflectional invariance along the x-axis. Expanding F in its arguments 
and keeping solely lowest order terms yields 

d t H = F + b(H -H) + a d x H + (a lx d 2 x + a ly d 2 y )H + 

+a 2 V 4 H + a 3x (d x H) 2 + a 3y (d y H) 2 + rj. (1) 

The physical significance of the terms on the rhs of (0) are constant erosion velocity Fq 
if all other terms on the rhs were zero, surface drift due to oblique incidence (ao-term), 



redeposition (6-term), anisotropic surface roughening via Bradley-Harper mechanism [h| 
(ai x , aiy terms), isotropic thermal surface diffusion (a2-term), anisotropic tilt-dependent 
sputter yield (a 3x , a 3y terms) and isotropic stochasticity in the erosion process with Gaussian 
white noise rj(x,y,t) of covariance 2D. 

Next we simplify (Q). First, the drift term a>od x H in is eliminated via the transforma- 
tion H (x, y, t) — > H(x — aot, y, t) implying a motion of the resulting pattern with a constant 
speed — cio in the x direction without any change of shape of the morphology of H. Second, 

is recast to a local field equation via the temporally nonlocal transformation [16] 

h = H -H+ (F /6)(l - e bt ) + e bt t d t H e - bt 'dt' (2) 

J o 

yielding 

d t h = bh + {a lx dl + a ly d 2 )h + a 2 V 4 h + a 3x (d x h) 2 + a 3y (d y h) 2 + 77. (3) 



Third, rescaling the time by — {a\ x / a 2 )t — > t, length scales by ya^/c^x — > x, height by 
— (a 3x /ai x )h — > h and noise by {a 2 a 3x / a\ x )rj — > r\ and introducing the three coefficients 
7 = {a 2 /a\ x )b > 0, a = a ly /a lx , (3 = a 3y /a 3x , leads to 

d t h = -(7 + d 2 x + a d 2 y + V 4 )/i + (d x h) 2 + /3(d y h) 2 + rj. (4) 

The anisotropic, damped KS equation (adKSE) (JH) constitutes our minimal model for ero- 
sion under normal and oblique incidence. The parameters a and (3 measure the relative 
anisotropies of the surface roughening and the sputter yield whereas the third parameter 

7 > the rescaled damping^ The limit 7 = is the aKSE 0, 0, whereas the limit 

a = l = /3 yields the idKSE 

Two properties of Eq.(J3J are of subsequent interest. (A) The parameter range a > 1 can be 
mapped to the range < a < 1 and vice versa by flipping the coordinate system, i.e. x — > y, 
y — > x and rescaling time by a 2 t — > t, length scales by ^ax — > x, height by (f3/a)h — > h and 
the noise by (f3/a 3 )r] 77 and transforming the coefficents, 7 = 7/a 2 , a = 1/a and (3 = 
This yields the functionally equivalent form d t h = — (7+c^+a d 2 +V i )h+(d x h) 2 +/3(d y h) 2 +ri. 
Consequently, the parameter sets (7, a, /3) with a < 1 and (7/a 2 , 1/a, with 1/a — a > 1 
lead to the same dynamics. (B) The unsealed version of Eq.(jlJ), i.e. Eq.tjHJ), possesses the 
invariance under the combined transformation h — > —h,a 3x —a 3x ,a 3y — > —a 3y . This 
implies that the corresponding morphologies are reflected about h — if the signs of a 3x and 
a 3a; are inverted, and, consequently, the role of mounds and valleys of the surface profile are 
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FIG. 1: Stability diagram of the flat front solution hpp = for n = 0. Markers indicate parameter 
values for subsequently used plots of the morphology. 

exchanged. In Eq. (|4*|) this invariance is hidden in the scaling. 

Since the coefficients in (0) depend on the angle of incidence and various other parameters 
such as ion flux and type, target material etc. and are not well known, we subsequently 
study how the properties of (j3J parametrically depend on the anisotropies and the damping 
without reference to a specific material and, by that, explore the pattern forming properties 
of Eq.gJ). 

Onset of pattern formation. - The non-stochastic limit of (J3J) obviously possesses a flat 
front solution hpp = /i(x, t) = corresponding to Hpp = F t as primary pattern. Its linear 
stability, determined by the growth or decay of small superimposed perturbations h(x, t) oc 
exp(ik-x+<Tt) in the linearized version of Q yields the growth rate cr(k) = —'-y+kl + aky — k 4 
as function of the wave-number k = (k x , k y ). The corresponding critical wave numbers are 
determined by (i) k x = and k y = Joi/2 or (ii) k x = \Jl/2 and k y = 0. Consequently, 
h(x, t) = is linearly unstable (i) with respect to plane wave perturbations in y-direction if 
7 < a 2 /4 and a > 1 or (ii) with respect to plane wave perturbations in x-direction if 7 < 1/4 
and a < 1. In Fig. ^the results of the linear stability analysis of hpp = are shown. Note 
the significant difference to the aKSE corresponding to the vertical line at 7 = where for 
a lx < and/or a\ y < no stable flat front solutions can be obtained. Similar arguments 
also apply to any nonlinear extension of the aKSE. 

Instability of perfect ripples. - Assuming rj = and a perfect ripple state h n (x,y,t) = 
h n (x,t) normal to the incidence direction (n-ripples) obeying d t h n = —(7 + d 2 + d^)h n + 
(d x h n ) 2 , the evolution of a slightly perturbed height function h(x,y,t) = h n (x,t) +p(y,t) in 
y-direction is governed by d t p = —{j+ady+d^)p+f3{d y p) 2 . Linearizing and inserting periodic 
perturbations p = e tk v +at yields the maximum growth rate o max = —7 + a 2 /4 implying that 
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FIG. 2: Snapshot of the surface morphologies for 7 = 0.24, (3 = 2.0 each and a = 1.0, 0.95, 0.8, 0.75 
(from left to right). Upper row: early stages at t = 200,620,400,400 (from left to right). Lower 
row: later stages at t = 10 4 each. 

n-ripple patterns are unstable against p-perturbation for 7 < a 2 /4 with a critical wave 
number k = \J a/2. Similarly, a perfect ripple state h t (x,y,t) = h t (y,t) tangential to the 
incidence direction (t-ripples) is unstable with respect to a slight perturbation p(x,t) in 
x-direction for 7 < 1/4 with a critical wave number k = yl/2. Consequently, such perfect 
t- or n-ripple patterns are always unstable in the aKSE. In the adKSE, however, they are, 
at least for these specific types of perturbations, stable in some parameter ranges. Since 
Ziberi et al. |7[ were able to experimentally produce basically perfectly straight ripples, the 
necessity of the damping term in our model equation seems to stringent. 

Transition from ripples to dots. - To numerically investigate how hexagonally arranged 
dot patterns can turn into ripples, we fix the anisotropy of the nonlinear term in J1J to a 
representative value (3 = 2 and the damping 7 to a value in the unstable regime of Iiff = 
but close to the instability, 7 = 0.24, and focus on a < 1. Numerical simulations have 
been performed using a finite difference method with periodic boundary conditions, a mesh 
size 200 x 200, spatial step size dx = 1, temporal step size dt = 0.01, an initially flat 
surface and an amplitude of the Gaussian white noise A = 5 x 10~ 4 corresponding to a noise 
covariance D = 4.16 x 10~ 6 . Varying the anisotropy a off from its isotropic limit a = 1 in 
simulations of (jlj), the following scenario as shown in Fig. El can be observed: (i) for a = 1 
the initially stochastically rough surface first evolves into a worm-like structures that turn 
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at later stages in a hexagonally arranged pattern with some defects, (ii) Lowering a to a 
value of 0.95, there is a drastic change in the preliminary stage of the evolution; slightly 
deformed n-ripples with some superimposed deformation structure in the y-direction occur. 
At later stages, however, the pattern stabilizes into a rather regular hexagonally arranged 
structure with even fewer defects than in (i). Consequently, even for a ^ 1 and (3^1 
hexagonally arranged dot patterns can be observed. Note, however, that these structures are 
not perfectly hexagonal since, as a result of any anistropy, the amplitudes of the hexagonal 
dots are no longer identical, (iii) Lowering a further to a value of 0.8, the preliminary stages 
are n-ripples with only minor modulations; surprisingly, the later stages consist of a highly 
regular arrangement of dots, by far more pronounced than in the case (i) or (ii). In our 
simulations, the evolution to the dot structure occurs via a successive periodic lacing up 
of the original n-ripples in the ^/-direction until a basically hexagonally arranged structure 
being slightly elongated in the y-direction is reached, (iv) Lowering a to a value of 0.75, 
the preliminary stages consisting of slightly modulated n-ripples stabilize for later stages to 
a basically perfect n-ripple state implying that a transition from long-time dot patterns to 
long-time ripple patterns occurs in the range 0.78 < a < 0.8. 

To investigate this transition quantitatively, we found that the temporal evolution of the 
surface roughness, w(t) = ((h — h) 2 ) 1 ^ 2 with the overbar and (...) denoting the spatial and 
ensemble average, shows that the cases (ii) and (iii) are quite distinct from the other two 
cases. This is because w(t) exhibits a double plateau behavior reflecting the transient ripple 
and the long-time dot structures with a turning point between these plateaus, as shown in 
the left panel of Fig. El Therefore, this turning point can be considered as an indicator for 
the evolution to long-time dot structures; its temporal position increases with lowering a 
until it vanishes at infinity. To substantiate this, we study the dependence of this turning 
point t c , depicted by lozenges in the right part of Fig. |3J as function of a. Also marked are 
the regions where preliminary states distinct from almost perfect ripples occur. As can be 
read off from this figure, these preliminary states generically develop into n-ripple structures 
followed subsequently by a transition to long-time dot structures if a < 0.93. For values of 
a smaller than about 0.78 only n-ripple structures are observable. This signals a bifurcation 
between long-time ripples and long-time dot structures at that value of a. To elucidate 
these results, this transition is not just restricted to the close vicinity of normal to oblique 
angle of incidence; it can also occur when the ratio a = a\ y ja\ x is accidentally close to 
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FIG. 3: Left panel: Temporal evolution of the surface roughness w. for 7 = 0.24, (3 = 2.0, 
a = 1.0,0.9,0.82,0.795,0.78 corresponding to the curves from top to bottom using an ensemble 
average over 40 runs. Right panel: Transition time t c to ripples and dots as function of a with 
fixed parameters 7 = 0.24 and (3 = 2.0. The regions where dots or ripples exist are delimitated by 
the turning points of the surface roughness w(t) as indicated in the left panel. For a > 0.93 only 
one turning point can be distinguished. 



unity for more oblique angles of incidence and, consequently, also hexagonally ordered dot 
structures show up in such parameter ranges. So far, the transition from ripples to dots, 
ckt has been obtained for a fixed value of 7. Numerical tests indicate that the transition 
line cvt(7) merges into the linear instability 7 = 0.25 as a — > 1. Finally, as a result of the 
afore- mentioned property (A) of an analogous scenario happens close the instability of 
hpp = for a > with the roles of n- and t-ripples being exchanged. 

Effect of damping. - To understand the role of the non-local term in (0) or the damping 
in (J3J), respectively, we show as a representative example in the first panel from the left 
of FigE] late stages of the surface morphology for the same values as in the third row of 
Fig. except that the damping parameter has been drastically reduced to a value of 
7 = 0.06. Obviously, lowering 7 has the effect that the rather periodic hexagonal dot 
structure seen for 7 = 0.24 has been coarsened into a pattern consisting of irregularly sized, 
droplet-like structures with a preferential elongation in y-direction. Flipping the sign of a 
for this parameter constellation, as shown in the second and third panel of Fig. HJ the initial 
and late stages of the surface evolution exhibit strong, irregular superimposed modulations 
of the n-ripples in the ^-direction. Consequently, the damping 7 acts as a order-disorder 
parameter, i.e. it gradually changes the pattern from the maximally disordered limit for 
7 = 0, i.e. the aKSE limit, to a quite regular ripple or dot pattern close to 7 = 0.25 (if 
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FIG. 4: Influence of the damping 7 on the surface morphology. From left to right: low damping 
7 = 0.06 for a late stage (t = 10 4 ) at a = 0.8, (3 = 2.0, an early stage (t = 50) at a = —0.8, 
(3 = 2.0, and corresponding late stage (t = 10 4 ); strong damping 7 = 0.3 for a late stage (t = 10 4 ) 
at a = 0.8, P = 2.0. Surface roughness for these four parameter constellations averaged over 40 
runs, left to right: w = 1.25, 1.02, 1.16, 1.79 • 10~ 3 . 

a < 1). For even larger 7, the deterministic limit of the adKSE possesses a linearly stable 
flat surface. Due to the nonzero stochasticity rj in our simulations, however, the resulting 
pattern for such values of 7 is stochastically rough, as depicted in the forth panel of Fig0] 
Most remarkably, the corresponding surface roughness is drastically reduced (by almost 
three orders of magnitude) in comparison to the cases with 7 < 0.24. Therefore, the pattern 
for 7 > 1/4, can be interpreted as a basically smooth surface with some stochastic variations 
that reflect the stochastic component of the sputter process. Such patterns have recently 
been reported by Frost et al. [17]. As a consequence, transitions to such basically smooth 
patterns triggered by changing the angle of incidence are preluded by a quite regular ripple 
or dot pattern. Assuming a nonlinear dependence of a and 7 on the angle of incidence and 
keeping (3 fixed (cf. also FigHJ), a variety of distinct transitions between patterns such as 
n-ripples <-> smooth surfaces <-> t- or n-ripples, n-ripples <-> dot structures <-» t- or n-ripples 
as well as more disordered states seem to be possible. Some of these transitions have been 



experimentally reported in (cf. also our introductory part). Using property (B) of 

Eq. (J3|). we also infer that patterns that look like negative images (e.g. mounds exchanged by 
holes) of the ones in Fig0]can be seen by going to negative values of a$ x and a^ y . Finally, our 
numerical simulations show that the transition to hexagonally arranged dots is generically 
subcritical implying that starting from an appropriately prestructured surface, dot patterns 
might evolve even for 7 > 1/4. 

Cancellation modes (CM). - Straight ripples oriented obliquely to the x-direction, as 
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found by Rost and Krug in the aKSE, can also be recovered in the adKSE. Assuming 
solutions h(x, y, t) = f(x—uy, t) that are constant along lines x = s+uy and demanding that 



the nonlinearity in the adKSE vanishes, implies u = ±y — 1/(3. Their evolution governed by 
dtf = — (7 + 1 + au 2 )d 2 f — (1 + u 2 ) 2 dgf has the maximum growth rate a max = —7 + {(3 — 



a) /[4(1 — (3) \ at a wave number q = J (3{(3 — a)/2({3 — l) 2 . Consequently, CM exist only 



if (3 < and a > (3 as in the case of the aKSE 



101 ]. Remarkably, CM exist as stationary 



solutions if 7 = (J3 - a) 2 /[A(l - (3) 2 ] in contrast to the aKSE 0. 
Conclusions. - The adKSE (0J) presented here seems to be a promissing candidate for a 
continuum model for low energy ion-beam erosion under oblique ion incidence since it can 
reproduce many distinct types of patterns also seen in experiments. To substantiate this, 
comparison with so far not available detailed experimental studies is needed. A more detailed 
discussion of the deterministic adKSE focussing on the bifurcation structure of its patterns 
will be given elsewhere |l8|. 
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